Analyse spatiale et territoriale
Formation Carthageo-Geoprisme 2023
Stage Carthageo-Geoprisme 2023
Objectifs de la semaine
- Développer les connaissances en analyse quantitative des données dans l’espace
- Exploiter des données individuelles (demandes de valeurs foncières, recensement de la population)
- Approfondir l’utilisation du logiciel R
- Travailler en groupe et préparer une restitution commune
- Répondre à une “commande” dans un temps limité
Organisation de la semaine
Planning des intervenants et contenu des journées
- 09/10 (matin, Pierre Pistre) : Organisation de la semaine & présentation des données
- 09/10 (après-midi, Claude Grasland) : Prise en main des données de recensement (fichier LOGEMT)
- 10/10 (matin, Pierre Pistre) : Prise en main des données de valeurs foncières (DVF)
- 10/10 (après-midi, Claude Grasland) : extractions des données pour la commande, et approfondissements statistiques & cartographiques
- 11/10 : Accompagnement des projets de groupe (Hadrien Commenges)
- 12/10 : Accompagnement des projets de groupe (Hadrien Commenges)
- 13/10 : Finalisation des projets de groupe (matin) et restitution orale (après-midi)
Horaires et planning des salles (Campus des Grands Moulins, Université Paris Cité)
- Horaires indicatifs : 9h30-12h30 et 13h30-17h.
- Salles : du 09/10 au 11/10 : salle de cours 209 (Bâtiment Olympe de Gouges, 2ième étage) ; 12/10 : salle informatique 375 (Bâtiment Olympe de Gouges, 3ième étage) ; 13/10 : salle de cours 209 (Bâtiment Olympe de Gouges, 2ième étage)
Données utilisées
(sources) Base de données principales :
- Fichier “Logements oridinaires” (LOGEMT) 2019 (localisations : IRIS, communes…), produit par l’INSEE à partir du recensement de la population. Données en accés libre : https://www.insee.fr/fr/statistiques/6544344?sommaire=6456104
- “Demandes de valeurs foncières” (version Etalab et extraction avec Opendatasoft), produit par le Ministère de l’Économie, des Finances et de la Souveraineté industrielle et numérique : https://www.data.gouv.fr/fr/datasets/demandes-de-valeurs-foncieres-geolocalisees/
(sources) Fichiers complémentaires :
- Table d’appartenance géographique des communes (INSEE) : https://www.insee.fr/fr/information/2028028
- Shapefile ADMIN-EXPRESS-COG® (IGN) des découpage administratifs de la France métropolitaine : https://geoservices.ign.fr/adminexpress#telechargement
- Shapefile Contours IRIS® (IGN) : https://geoservices.ign.fr/contoursiris
Exercice d’application en groupe
Périmètres de l’exercice
- Objet d’étude : les dynamiques immobilières dans les principales agglomérations administratives de la “mégarégion parisienne” (hors Paris)
- Population d’étude : les logements vendus récemment et leurs habitants
- Espace(s) d’étude : les principales agglomérations ayant notamment le statut de “métrople” ou de “communauté urbaine” (hors Paris)
- Echelles et mailles géographiques : département, commune, IRIS
Consignes de la “commande”
- Contexte général : il est connu que Paris influence la dynamique des territoires bien au-delà de ses territoires les plus proches, à commencer par les principales agglomérations urbaines de la “mégarégion parisienne” (au sens élargi de : https://atlas-paris-mega-region.univ-rouen.fr/). Le secteur immobilier est particulièrement concerné du fait d’une progression sensible des prix des logements dans la plupart de ces territoires au cours à minima de la dernière décennie et d’un potentiel renforcement avec la crise COVID par une progression des installations résidentielles. Pour autant, la géographie des dynamiques immobilières internes à ces agglomérations urbaines reste mal connue de même que les facteurs explicatifs (influence parisienne, spécificités locales et régionales…) et les tendances récentes.
- Commande : après avoir fait le constat de dynamiques immobilières souvent proches, les principales agglomérations urbaines de la “mégarégion parisienne” ont décidé de s’associer pour demander aux différentes antennes régionales de l’Insee la réalisation d’études rigoureuses et fines spatialement sur la situation immobilière de leurs différents territoires de compétence et périphéries proches. Le bilan de ces études sera présenté dans des publications synthétiques, intelligibles par le plus grand nombre et par les acteurs concernés, sous la forme de plusieurs 4 pages - comme produits régulièrement par l’Insee pour publier les résultats de ses études (par exemple, file:///home/pierre/T%C3%A9l%C3%A9chargements/ip1715.pdf ; file:///home/pierre/T%C3%A9l%C3%A9chargements/IR134_Notaires-IPLA_1T-23.pdf).
Modalités de l’exercice, restitution du travail et organisation durant la semaine
- Constitution de 6 (voire 7) groupes de 4-5 étudiant.e.s, mélangeant les profils de Master (= Lundi matin)
- Réflexion sur un angle problématique pour chaque cas d’étude, ainsi que les possibilités d’analyse et d’exploration des données : variables pertinentes, traitements envisagés… (= Lundi à Mercredi)
- Réalisation des analyses statistiques et cartographiques (= Mardi à Jeudi)
- Organisation, mise en forme des analyses et rédaction de la note écrite (= Mardi à Vendredi matin)
- Présentation orale (environ 10 minutes, sans autre support visuel que la note écrite) devant un “jury” composé des intervenants de la semaine (= Vendredi après-midi)
Cas d’étude
- Caen (communauté urbaine)
- Le Mans (communauté urbaine)
- Orléans (métropole)
- Reims (communauté urbaine)
- Rouen (métropole)
- Tours (métropole)
En option supplémentaire : Amiens (communauté d’agglomération)
Format du rendu (type 4 pages Insee)
- Introduction et définition (thématique et statistique) de l’objet : thème et cas d’étude
- Etat des lieux général des prix immobiliers dans le cas d’étude (données : DVF ; échelle : “agglomération” (rayon 40km) ; maille : commune)
- Approfondissements sur les prix immobiliers par type de biens ou de localisation (données : DVF ; échelles : “agglomération” (rayon 40km) et zooms ; mailles : commune, IRIS)
- Profils des acheteurs au sens des nouveaux arrivants (moins de 3 ans) qui sont propriétaires (données : recensement de la population ; échelles : “agglomération” (rayon 40km) et zooms ; mailles : commune, IRIS)
RP2 : Agrégation territoriale
Le format sf (spatial features)
La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois
- un tableau de données (l’équivalent du fichier .dbf)
- une géométrie (l’équivalent du fichier .shp)
- une projection (l’équivalent du fichier .prj)
Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou
dans d’autres formats standards comme GeoJson, la première tâche
consiste donc à les convertir au formt sf afin de pouvoir les utiliser
facilement dans R. L’importation se fait à l’aide de l’instruction
st_read en indiquant juste le nom du fichier .shp à
charger. Les autres fichiers (.dbf ou .proj) seront lus également et
intégrés dans l’objet qui hérite de la double classe data.frame
et sf.
Etapes de préparation des données
Dans notre exemple, nous allons suivre les étapes suivantes :
- Préparer les données statistiques par IRIS dans un data.frame
- Charger un fonds de carte par IRIS au format sf
- Effectuer une jointure entre les deux fichiers par le code IRIS
- Sauvegarder le résultat
- Agréger les données statistiques et géométriques par commune
- Sauvegarder le résultat.
Préparer les données statistiques
On importe le fichier des individus et on corrige le code IRIS des communes non divisées en IRIS.
myrep<-"data/Rennes/"
tab_ind<-readRDS(paste0(myrep,"EPCI_men.RDS"))
tab_ind$IRIS[tab_ind$IRIS=="ZZZZZZZZZ"]<-
paste0(tab_ind$INSEE_COM[tab_ind$IRIS=="ZZZZZZZZZ"],"0000")
head(tab_ind[,1:5],3)FALSE COMMUNE ARM IRIS ACHL AEMM
FALSE 1: 35001 ZZZZZ 350010101 B12 1976
FALSE 2: 35001 ZZZZZ 350010101 B12 1977
FALSE 3: 35001 ZZZZZ 350010101 B12 0
Selectionner les lignes et colonnes
On décide de ne sélectionner que les ménages de propriétaires installés depuis moins de 3 ans et on retient comme variable le type de logement que l’on recode en deux catégories
tab_ind2 <- tab_ind %>% filter(as.numeric(ANEM) < 3, # durée d'occupation < 3 ans
STOCD =="10", # propriétaire de son logement
TYPL %in% c(1,2)) %>% # Maison ou Appartement
mutate(TYPL = as.factor(TYPL)) %>%
select(IPONDL, IRIS, TYPL)
levels(tab_ind2$TYPL)<-c("Maison","Appt")
knitr::kable(head(tab_ind2,5),digits=2)| IPONDL | IRIS | TYPL |
|---|---|---|
| 1.07 | 350010101 | Maison |
| 1.03 | 350010101 | Maison |
| 1.07 | 350010101 | Maison |
| 0.99 | 350010101 | Maison |
| 1.11 | 350010101 | Maison |
Agréger les données
On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :
tab_long<- tab_ind2 %>%
group_by(IRIS,TYPL)%>%
summarise(NB=sum(IPONDL))
knitr::kable(head(tab_long,5),digits=2)| IRIS | TYPL | NB |
|---|---|---|
| 350010101 | Maison | 36.13 |
| 350010101 | Appt | 9.71 |
| 350010102 | Maison | 106.53 |
| 350010102 | Appt | 57.60 |
| 350220000 | Maison | 21.89 |
Pivoter le tableau
Puis on fait “pivoter” le tableau pour l’obtenir en format large :
tab_large <- tab_long %>% pivot_wider(id_cols = IRIS,
names_from = TYPL,
names_prefix = "T_",
values_from = NB,
values_fill = 0)
knitr::kable(head(tab_large,5),digits=2)| IRIS | T_Maison | T_Appt |
|---|---|---|
| 350010101 | 36.13 | 9.71 |
| 350010102 | 106.53 | 57.60 |
| 350220000 | 21.89 | 0.99 |
| 350240101 | 142.85 | 41.34 |
| 350240102 | 34.58 | 10.50 |
Ajouter de nouvelles variables
On ajoute de nouvelles variables telles que le nombre total de nouveaux propriétaires et la part des propriétaires de maisons parmi eux.
tab<- tab_large %>% mutate(T_TOT = T_Maison+T_Appt,
Maison_pct = 100*T_Maison/T_TOT)
knitr::kable(head(tab,5),digits=c(0,0,0,0,1))| IRIS | T_Maison | T_Appt | T_TOT | Maison_pct |
|---|---|---|---|---|
| 350010101 | 36 | 10 | 46 | 78.8 |
| 350010102 | 107 | 58 | 164 | 64.9 |
| 350220000 | 22 | 1 | 23 | 95.7 |
| 350240101 | 143 | 41 | 184 | 77.6 |
| 350240102 | 35 | 10 | 45 | 76.7 |
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de nouveux propriétaires en maison individuelle.
programme
p <- ggplot(tab) + aes (x = Maison_pct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90, 100)) +
scale_x_continuous("% de Maisons") +
scale_y_continuous("Nombre d'IRIS") +
ggtitle(label = "Type de logement des nouveaux propriétaires",
subtitle = "Source : INSEE, RP 2019 - Métropole de Rennes")
pCharger les données géométriques
On importe le fichier des iris de l’agglomération qui est de type sf (spatial features)
map_iris <- readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
map_iris<-map_iris[,c(5,6,3,2)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")
class(map_iris)FALSE [1] "sf" "data.frame"
| IRIS | NOM_IRIS | COM | NOM_COM | geometry |
|---|---|---|---|---|
| 352381204 | Les Cloteaux | Rennes | 35238 | MULTIPOLYGON (((350568 6786… |
| 352750000 | Saint-Gilles | Saint-Gilles | 35275 | MULTIPOLYGON (((337255.9 67… |
Visualisation du fonds iris avec sf
On peut facilement faire un plot de la colonne geometry du fichier sf
Jointure des données IRIS et du fonds de carte
programme
map_iris_tab<-merge(map_iris,tab,
by.x="IRIS",by.y="IRIS",
all.x=T,all.y=F)
knitr::kable(head(map_iris_tab,3),digits=2)| IRIS | NOM_IRIS | COM | NOM_COM | T_Maison | T_Appt | T_TOT | Maison_pct | geometry |
|---|---|---|---|---|---|---|---|---|
| 350010101 | Ouest | Acigné | 35001 | 36.13 | 9.71 | 45.84 | 78.82 | MULTIPOLYGON (((360989.4 67… |
| 350010102 | Est | Acigné | 35001 | 106.53 | 57.60 | 164.13 | 64.90 | MULTIPOLYGON (((364373.3 67… |
| 350220000 | Bécherel | Bécherel | 35022 | 21.89 | 0.99 | 22.88 | 95.65 | MULTIPOLYGON (((333530.4 68… |
Cartographie rapide avec sf
Une astuce pour visualiser rapidement une carte choropléthe avec le seul packahe sf :
Sauvegarde du fichier par IRIS
Agrégation statistique + géométriques
Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”
Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.
Agrégation des IRIS en communes
L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries
Agrégation des iris en communes
résultat statistique
| COM | NOM_COM | T_Maison | T_Appt | T_TOT | Maison_pct |
|---|---|---|---|---|---|
| Acigné | 35001 | 143 | 67 | 210 | 67.9 |
| Betton | 35024 | 291 | 82 | 373 | 78.0 |
| Bourgbarré | 35032 | 162 | 24 | 186 | 87.2 |
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de nouveaux propriétaires en maison par commune.
p <- ggplot(map_com_tab) +aes (x = Maison_pct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90, 100)) +
scale_x_continuous("% de Maisons") +
scale_y_continuous("Nombre de communes") +
ggtitle(label = "Type de logement des nouveaux propriétaires",
subtitle = "Source : INSEE, RP 2019, Agglomération de Rennes")
p RP3 : Cartographie statique
Le package map_sf
Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.
On trouvera la documentation du package mapsf à l’adresse suivante :
Création d’un template cartographique
Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :
tracé d’un fonds de carte vierge
La fonction mf_map() avec le paramètre
type = "base"permet de tracer une carte vide
Superposition de couches
On peut toutefois ajouter toute une série de paramètres
supplémentaire (col=, border=,
lwd=, …) et superposer plusieurs fonds de carte avec le
paramètre add = TRUE. L’ajout de la fonction
layout permet de rajouter un cadre une légende.
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray50",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=0.6,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude",
credits = "Sources : IGN et INSEE")Ajout d’un thème
On peut finalement modifier l’ensemble de la carte en lui ajoutant
une instruction mf_theme() qui peut reprendre des styles
existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”,
“candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”,
“barcelona”) mais aussi créer ses propres thèmes
#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray50",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=0.6,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude",
credits = "Sources : IGN et INSEE")Ajout de texte
On peut ajouter une couche de texte avec la fonction
mf_label(). Par exemple, on va ajouter à la carte
précédente le nom des communes
programme
mf_theme("agolalight")
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray50",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=0.6,
add = TRUE)
# Ajoute une couche de labels
mf_label(map_com,
var="COM",
cex=0.4,
col="blue",
overlap = FALSE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude",
credits = "Sources : IGN et INSEE")Carte de stock
Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.
Dans le package mapsf, on réalise ce type de carte à
l’aide de la fonction mf_map()en lui donnant le paramètre
type="prop".
On va tenter à titre d’exemple de représenter la distribution du nombre de nouveaux propriétaires.
Carte de stock minimale
Les instructions minimales sont les suivantes :
# Trace les contours des communes
mf_map(x= map_iris,
type = "base")
# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris,
type ="prop",
var = "T_TOT",
add=TRUE)Mais le résultat est peu satisfaisant car les cercles sont trop
grands. Il faut en pratique toujours effectuer un réglage de ceux-ci
avec l’instruction inches=
Carte de stock habillée
mf_theme("agolalight")
mf_map(map_iris, type = "base",
col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="black",lwd=1,add = TRUE)
mf_map(map_iris, var = "T_TOT",type = "prop",
inches = 0.05, col = "red",leg_pos = "left",
leg_title = "Nombre de ménages", add=TRUE)
mf_layout(title = "Distribution des nouveaux propriétaires",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte choroplèthe
Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.
La fonction du package mapsf adaptée aux variables
d’intensité est la fonction mf_map()munie du paramètre
type = "choro".
On va prendre l’exemple du % de maisons parmi les nouveaux propriétaires
Carte choroplèthe minimale
Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).
Carte choroplèthe habillée
On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.
mybreaks = c(0, 10,20,30,40,50,
60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5),
pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris, var = "Maison_pct",
type = "choro",
breaks = mybreaks,pal = mypal,
border=NA,
col_na = "gray80",leg_title = "% maisons",
leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com, type = "base", col = NA,
border="black",lwd=1,add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte stock + choroplèthe (1)
On peut combiner les deux modes cartographiques par superposition :
mf_theme("agolalight")
# Choisit les classes
mybreaks = c(0,20,40,60,80,100)
# Trace la carte choroplèthe
mf_map(
x = map_iris,
var = "Maison_pct",
breaks = mybreaks,
# pal=mypal,
type = "choro",
border="white",
col_na = "gray80",
lwd=0.3,
leg_title = "% Maisons",
leg_val_rnd = 0,
)
# Ajoute les cercles proportionnels
mf_map(
x =map_iris,
var = "T_Maison",
type = "prop",
inches = 0.06,
col = "red",
leg_pos = "right",
leg_title = "Effectif",
add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="black",
lwd=1,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les nouveaux propriétaires de maison",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte stock + choroplèthe (2)
Mais les cercles dissimuent alors les plages de couleur, aussi on
peut utiliser le type prop_choro qui place la variable
choroplèthe à l’intérieur des cercles
mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens","Reds"))
mf_map(map_iris, type = "base",
col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris, var = c("T_TOT", "Maison_pct"),
inches = 0.06, col_na = "grey", pal=mypal,
breaks = mybreaks, nbreaks = 4, lwd = 0.1,
leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
leg_title = c("Nouveaux propriétaires", "% Maison"),
add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons",
frame = TRUE, credits = "Sources : IGN et INSEE")DVF1 : Acquisition des données
Introduction
Objectif
On se propose de collecter l’ensemble des ventes de maison ou d’appartement dans un rayon de 50 km autour d’une ville. On choisit ici Rennes en prenant la place du Parlement de Bretagne comme centre de référence
Données dvf
Les données sur les dvf sont disponibles sur le site publicopendatasof:
https://public.opendatasoft.com
On s’intéresse plus particulièrement aux données dvf géolocalisées accessibles ici :
Sélection des informations
Le fichier comporte 28 millions d’enregistrement ce qui est évidemment beaucoup … On va donc procéder à une sélection en se servant des différents onglets disponibles.
A titre d’exemple, nous allons essayer de télécharger les enregistrements vérifiant les conditions suivantes :
- localisation dans un rayon de 60 km autour d’Amiens
- type de transaction : Vente
- type de bien : Maison ou appartement
Nous constatons qu’il y a 133294 enregistrements vérifiant ces conditions.
Lien de téléchargement
Nous pouvons récupérer les données soit au format .csv
pour les analyse statistiques, soit au format .geojson pour
les analyses spatiales. Mais dans ce cas on ne peut pas dépasser la
taille de 50 000 enregistrement. Il vaut donc mieux télécharger au
format .csv puisque ce fichier comporte les latitudes et longitudes ce
qui nous permettra de créer un fichier cartographique a posteriori
On peut effectuer le téléchargement depuis le site web ou bien juste enregistrer le lien de téléchargement et effectuer l’opération dans R.
Pour trouver le lien on passe la souris au dessus du lien “télécharger les données au format .csv / seulement les enregistrements sélectionnées” et on effectue un click droit suivi de “copier le lien”
Construction d’une API
On récupère le lien de téléchargement du fichier .csv :
A première vue ceci paraît assez obscur mais on réalise assez vite que cette URL reprend l’ensemble des spécifications envoyées pour sélectionner nos enregistrements. On peut souligner en gras les parties utiles :
- &refine=nature_mutation%3A%22Vente%22&
- &refine=type_local%3A%22Maison%22
- &refine=type_local%3A%22Appartement%22
- &where=(distance(%60geo_point%60%2C%20geom%27
- POINT(2.26968%2049.90807)%27%2C%2060037m%20))
Nous allons donc créer un lien qui extrait automatiquement les dvf en modifiant juste la position du point central et le rayon. On prend cette fois-ci nos paramètres de Rennes
myurl<-paste0("https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/exports/csv?lang=fr&refine=nature_mutation%3A%22Vente%22&refine=type_local%3A%22Maison%22&refine=type_local%3A%22Appartement%22&facet=facet(name%3D%22nature_mutation%22%2C%20disjunctive%3Dtrue)&facet=facet(name%3D%22type_local%22%2C%20disjunctive%3Dtrue)&where=(distance(%60geo_point%60%2C%20geom%27POINT(", latctr,"%20",lonctr,")%27%2C%20",rayon,"m))&timezone=Europe%2FParis&use_labels=true&delimiter=%3B")Téléchargement
On lance la requête avec l’instruction download.file()
et on stocke le résultat dans notre répertoire de travail :
essai de l’URL ‘https://public.opendatasoft.com/…’ downloaded 38.4 MB
L’opération s’execute en un peu moins d’une minute dans cet exemple. Mais il peut arriver qu’elle échoue lorsque le serveur est saturé et dans ce cas il faut recommencer.
DVF2 : Nettoyage des données
Introduction
Objectif
On se propose de nettoyer les donénes relative à l’ensemble des ventes de maison ou d’appartement dans un rayon de 50 km autour d’une ville. On rappelle les paramètres qui ont été utilisés pour acquérir ces données.
Paramètres
Merci à Boris et Florent !
Nous allons effectuer le nettoyage en nous appuyant sur les propositions de Boris Mericskay et Florent Demoraes parues dans https://journals.openedition.org/cybergeo/39583
Les auteurs fournissent en effet un accès à tous les programmes utilisés dans leur article via une archive github : https://github.com/ESO-Rennes/Analyse-Donnees-DVF
Chargement du jeu de données
On adapte le programme de Boris et Florent en utilisant la fonction
fread() du package data.table car elle est
plus rapide/
Etape 1 : Sélection des mutations de type “Ventes” de “Maison” et “Appartement’
Programme
Remarque
Etape inutile si on a déjà fait ce double filtrage à l’aide de l’API
Etape 2 : Sélection et renommage des variables
etape2 <- etape1bis %>%
select(id = `Identifiant de mutation (Etalab)`,
disposition = `Numéro de disposition` ,
parcelle= `Identifiant de la parcelle cadastrale`,
date = `Date de la mutation`,
nature = `Nature de la mutation`,
INSEE_COM = `Code INSEE de la commune`,
INSEE_DEP = `Code INSEE du département` ,
type = `Type de local` ,
surface = `Surface réelle du bâti` ,
piece = `Nombre de pièces principales`,
prix = `Valeur foncière` ,
latitude = `Latitude`,
longitude = `Longitude`)Etape 3 : Remplacement des cellules vides par des NA et suppression des NA
Etape 4 : Sélections des mutations simples
Regrouper les transactions selon l’ID, la surface et la vente
Etape 5 : Jointure attributaire pour récupérer les informations de la mutation
jointure
Modification des formats des colonnes
Suppression des valeurs aberrantes
Fixer un seuil minimal et maximal des prix (percentiles)
1% 99%
17000 616320
Fixer un seuil minimal et maximal des prix (histogramme)
Créer deux jeux de données (maisons / appartements)
Appartement Maison
42003 68734
Etape 6 : Sélection des bornes de prix et de surface
Etape 8 : Sélection des bornes de prix au m2
Fixer un seuil minimal des prix au m2 (percentile)
1% 99%
387.3057 5825.8048
Etape 9 : Corrections diverses
Transformation de la date en année
Arrondir les variables numériques
Etape 10 : Exportation
Mise en ordre des lignes et colonnes
Ecrire le jeu de données final en csv
DFVF 3 : Cartographie
Retour sur sf
Nous revenons sur le package sf (spatial features)
que nous avons déjà rencontré au moment de la création de cartes
thématiques par IRIS ou communes à l’aide du package
mapsf.
Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des ventes foncières au dessus du fonds de carte des IRIS ou communes.
Mais il pourra aussi servir de base à des cartographies
dynamiques permettant de placer les points correspondant aux
ventes sur des réseaux de rue et plus généralement sur des “tuiles”
cartographiques permettant d’effectur des zoom. On utilisera à cet effet
d’autres packages comme leaflet ou sa version simplifiée
mapview.
Zone d’étude
On définit les paramètres de la zone d’étude en ajoutant le lien vers les dvf nettoyées
Choix de communes
On choisit plus précisément de travailler sur les communes de l’agglmomération
map_com <- readRDS(paste0(myrep,"EPCI_mapcom.RDS")) %>%
select(INSEE_COM, NOM, geometry)
plot(map_com$geometry,
col="lightyellow")Sélection des dvf nettoyées
Nous reprenons le fichier nettoyé de localisation des ventes de maison et d’appartement établi au chapitre précédent et nous ne conservons que 9 variables ainsi que les ventes situées dans les communes de la zone choisie
dvf <- readRDS(paste0(myrep,mydvf)) %>%
select(annee, INSEE_COM,
type, prix, surface, prixm2,
latitude, longitude) %>%
filter(INSEE_COM %in% map_com$INSEE_COM)
kable(head(dvf,3))| annee | INSEE_COM | type | prix | surface | prixm2 | latitude | longitude |
|---|---|---|---|---|---|---|---|
| 2014 | 35238 | Appartement | 81000 | 66 | 1227 | 48.11526 | -1.692243 |
| 2014 | 35238 | Appartement | 98500 | 33 | 2985 | 48.11689 | -1.681832 |
| 2014 | 35238 | Appartement | 121000 | 62 | 1952 | 48.10974 | -1.705489 |
Transformation des dvf en fichier sf
On ajoute une colonne geometry en utilisant la fonction
st_as_sf()du package sf. On attribue à la
projection le code EPSG = 4326 qui correspond à un fichier latitude
longitude.
map_dvf <- st_as_sf(dvf, coords = c("longitude","latitude"))
st_crs(map_dvf)<- 4326
kable(head(map_dvf,3))| annee | INSEE_COM | type | prix | surface | prixm2 | geometry |
|---|---|---|---|---|---|---|
| 2014 | 35238 | Appartement | 81000 | 66 | 1227 | POINT (-1.692243 48.11526) |
| 2014 | 35238 | Appartement | 98500 | 33 | 2985 | POINT (-1.681832 48.11689) |
| 2014 | 35238 | Appartement | 121000 | 62 | 1952 | POINT (-1.705489 48.10974) |
Test de superposition
On tente de superposer les deux cartes des ventes et des communes.
Commune puis dvf
Vérification de la projection
Nous savons que les coordonnées du fichier DVF sont en latitude longitude (EPSG=4326). Mais la projection de notre fonds des communes est en Lambert-93 (EPSG = 2154)
Coordinate Reference System:
User input: RGF93 v1 / Lambert-93
wkt:
PROJCRS["RGF93 v1 / Lambert-93",
BASEGEOGCRS["RGF93 v1",
DATUM["Reseau Geodesique Francais 1993 v1",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4171]],
CONVERSION["Lambert-93",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",46.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",3,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",49,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",44,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",700000,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",6600000,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["France - onshore and offshore, mainland and Corsica."],
BBOX[41.15,-9.86,51.56,10.38]],
ID["EPSG",2154]]
A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS
Ajustement des deux projections
On transforme donc le fichier des DVF en Lambert-93 pour obtenir une adéquation avec le fichier des communes.
Et cette fois-ci on réussit à superposer les deux fonds de carte
Cartographie ponctuelle avec mapsf
On peut utiliser les fonctions de mapsf vues précédemment. Par exemple repérer les ventes de maison et d’appartement :
mf_map(map_com, type="base",
col="lightyellow")
mf_map(map_dvf, type = "typo",
var= "type",add =T,
border=NA,
cex=0.3,
leg_title = "Bien vendu")
mf_label(map_com,
var="NOM",cex = 0.4,
overlap = F,
col="black",
halo=T)
mf_layout(title = "Ventes de maisons et d'appartement 2013-2020",
credits = "Sources : IGN, INSEE, DVF",
scale = T)Cartographie ponctuelle avec mapview
Le package mapview pemet de créer facilement des cartes interactives ù l’on pourra visualiser les ventes sur un fonds de carte de son choix (OSRM, ESRI, …). On peut l’utiliser par exemple pour voirune commune précise
On peut améliorer sensiblement la visualisation à l’aide de paramètres présentés sur le site web du package https://r-spatial.github.io/mapview/
library(mapview)
# Choix de la zone d'étude
map_com_sel<-map_com %>% filter(INSEE_COM == "35196")
map_dvf_sel<-map_dvf %>% filter(INSEE_COM == "35196")
# Choix des tuiles
mapviewOptions(basemaps = c("OpenStreetMap" ,
"Esri.WorldImagery"))
# Première couche
map1 <- mapview(map_com_sel,
col.regions="lightyellow",
alpha.regions = 0.4)
# Deuxième couche
map2<-mapview(map_dvf_sel,
zcol = "type",
cex = 4,
col.regions=c("blue","red"),
alpha.regions = 0.4)
# Assemblage
map1+map2Agrégation par communes
On peut évidemment agréger les informations sur les valeurs foncières par communes. Par exemple compter le nombre de maisons vendues et leur prix moyen au m2
dvf_by_com <- dvf %>%
group_by(INSEE_COM) %>%
filter(type == "Maison") %>%
summarise(nb = n(),
med_prixm2 = median(prixm2))
kable(head(dvf_by_com,3))| INSEE_COM | nb | med_prixm2 |
|---|---|---|
| 35001 | 400 | 2417.0 |
| 35022 | 66 | 1369.5 |
| 35024 | 745 | 2576.0 |
On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :
On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :
mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_com_dvf, type = "choro",
var = "med_prixm2",
breaks=mybreaks, pal = mypal,
border="gray", lwd=0.3,
leg_title = "prix en €/m2", leg_val_rnd = 0,
leg_pos = "left")
mf_map( map_com_dvf, type = "prop",
var = "nb", inches = 0.12,
col="black",border="white",
leg_title = "nb. de ventes",
leg_pos = "topleft")
mf_layout(title = "Les ventes de maison de 2013 à 2020",
frame = TRUE, credits = "Sources : IGN et DVF",
arrow = F )Cartographie par IRIS
Supposons maintenant que l’on souhaite cartographier les prix des maisons par IRIS à l’intérieur de la commune de Rennes. Comment faire dans la mesure où l’IRIS n’est pas mentionné dans le fichier DVF ? Réponse : en superposant les deux fonds de carte et en les intersectant.
Superposition
map_iris_ctr<-readRDS(paste0(myrep,"CTR_mapiris.RDS")) %>%
select("NOM_COM","INSEE_COM","CODE_IRIS","NOM_IRIS","geometry")
map_dvf_ctr <- map_dvf %>% filter(INSEE_COM %in% map_iris_ctr$INSEE_COM)
plot(map_iris_ctr$geometry, col="lightyellow")
plot(map_dvf_ctr$geometry,col="red",pch=21,cex=0.1,add=T)Intersection
On peut faire l’opération avec un SIG ou bien avec le package
sf en utilisant sa fonction st_intersect()
qui va ajouter tous les attributs du fichier IRIS au fichier DVF :
| annee | INSEE_COM | type | prix | surface | prixm2 | NOM_COM | INSEE_COM.1 | CODE_IRIS | NOM_IRIS | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 62 | 2014 | 35238 | Appartement | 80000 | 44 | 1818 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350454.3 6786503) |
| 130 | 2014 | 35238 | Appartement | 70000 | 55 | 1273 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350613.6 6786410) |
| 158 | 2014 | 35238 | Appartement | 115000 | 79 | 1456 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350694.4 6786355) |
| 212 | 2014 | 35238 | Appartement | 109000 | 73 | 1493 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350662.4 6786511) |
Agrégation
Il ne reste plus qu’à refaire l’agrégation en suivant la même procédure que pour les communes mais en retirant les iris où il y a eu moins de 5 ventes.
map_iris_dvf <- inter %>%
st_drop_geometry() %>%
group_by(CODE_IRIS, NOM_IRIS) %>%
filter(type == "Maison") %>%
summarise(nb = n(),
med_prixm2 = median(prixm2)) %>%
filter(nb >=5) %>%
right_join(map_iris_ctr) %>%
st_as_sf()
kable(head(map_iris_dvf,3))| CODE_IRIS | NOM_IRIS | nb | med_prixm2 | NOM_COM | INSEE_COM | geometry |
|---|---|---|---|---|---|---|
| 352380101 | Cathédrale | 5 | 2680 | Rennes | 35238 | MULTIPOLYGON (((351427 6789… |
| 352380102 | Hoche | 15 | 4671 | Rennes | 35238 | MULTIPOLYGON (((352032.6 67… |
| 352380109 | Vieux Saint-Étienne | 17 | 5263 | Rennes | 35238 | MULTIPOLYGON (((351501.3 67… |
Cartographie
Il ne reste plus qu’à faire la carte en reprenant le programme utilisé pour les communes :
mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_iris_dvf, type = "choro",
var = "med_prixm2",
breaks=mybreaks, pal = mypal,
border="gray", lwd=0.3,
leg_title = "prix en €/m2", leg_val_rnd = 0,
leg_pos = "left")
mf_map(map_iris_dvf, type = "prop",
var = "nb", inches = 0.08,
col="black",border="white",
leg_title = "nb. de ventes",
leg_pos = "topleft")
mf_layout(title = "Les ventes de maison de 2013 à 2020",
frame = TRUE, credits = "Sources : IGN et DVF",
arrow = F )